/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "fastonebigheader.h"

static inline float
squareit (float x)
{
  return x * x;
}

/* Ugh ... I should have used C++ */

#define logprating(_rating, _loga, _rank,                               \
                   _logbeta, _dlogbeta, _ddlogbeta, _c, _d, _n_labels)  \
  ({                                                                    \
    /*                                                                  
     * ok we want to compute the logarithm of
     * 
     * exp (exponents[rating])     exp (exponents[rating] - max)
     * ----------------------- = --------------------------------
     * sum_i exp (exponents[i])   sum_i exp (exponents[i] - max)
     *
     * the RHS is more numerically stable
     */                                                                 \
                                                                        \
    float _a[_rank];                                                    \
    for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)          \
      {                                                                 \
        _a[_logprn] = fastexp (_loga (_logprn));                        \
      }                                                                 \
    float _beta = fastexp (_logbeta);                                   \
    float _logp[_n_labels];                                             \
    unsigned int _maxindex = 0;                                         \
    float _invrank = 1.0 / _rank;                                       \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        _logp[_logpri] = 0;                                             \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _logp[_logpri] -= _a[_logprn] *                             \
                              _c (_logprn) *                            \
                              _d (_logprn, _logpri);                    \
          }                                                             \
                                                                        \
        _logp[_logpri] *= _beta * _invrank;                             \
                                                                        \
        if (_logp[_logpri] > _logp[_maxindex])                          \
          {                                                             \
            _maxindex = _logpri;                                        \
          }                                                             \
      }                                                                 \
                                                                        \
    float _max = _logp[_maxindex];                                      \
    float _denominator = 0;                                             \
    float _ddenominator_dlogbeta = 0;                                   \
    float _dddenominator_ddlogbeta = 0;                                 \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        float d = _logp[_logpri] - _max;                                \
        float expd = fastexp (d);                                       \
        _denominator += expd;                                           \
        _ddenominator_dlogbeta += expd * d;                             \
        _dddenominator_ddlogbeta += expd * d * (d + 1.0);               \
      }                                                                 \
                                                                        \
    _dlogbeta = _logp[_rating] - _max                                   \
               - _ddenominator_dlogbeta / _denominator;                 \
                                                                        \
    _ddlogbeta = _logp[_rating] - _max                                  \
                + squareit (_ddenominator_dlogbeta / _denominator)      \
                - _dddenominator_ddlogbeta / _denominator;              \
                                                                        \
    _logp[_rating] - _max - fastlog (_denominator);                     \
  })

#define dlogpratingdall(_pzj, _dloga, _dc, _dd, _rating,                \
                        _loga, _rank, _logbeta, _c, _d, _n_labels)      \
  ({                                                                    \
    float _a[_rank];                                                    \
    for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)          \
      {                                                                 \
        _a[_logprn] = fastexp (_loga (_logprn));                        \
      }                                                                 \
                                                                        \
    float _beta = fastexp (_logbeta);                                   \
    float _logp[_n_labels];                                             \
    unsigned int _maxindex = 0;                                         \
    float _invrank = 1.0 / _rank;                                       \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        _logp[_logpri] = 0;                                             \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _logp[_logpri] -= _a[_logprn] *                             \
                              _c (_logprn) *                            \
                              _d (_logprn, _logpri);                    \
          }                                                             \
                                                                        \
        _logp[_logpri] *= _beta * _invrank;                             \
                                                                        \
        if (_logp[_logpri] > _logp[_maxindex])                          \
          {                                                             \
            _maxindex = _logpri;                                        \
          }                                                             \
      }                                                                 \
                                                                        \
    float _max = _logp[_maxindex];                                      \
    float _denominator = 0;                                             \
    float _ddenominator_dloga[_rank][_n_labels];                        \
    float _ddenominator_dlogbeta = 0;                                   \
    float _ddenominator_dc[_rank][_n_labels];                           \
    float _ddenominator_dd[_rank][_n_labels];                           \
                                                                        \
    memset (_ddenominator_dloga,                                        \
            0,                                                          \
            sizeof (_ddenominator_dloga));                              \
                                                                        \
    memset (_ddenominator_dc,                                           \
            0,                                                          \
            sizeof (_ddenominator_dc));                                 \
                                                                        \
    memset (_ddenominator_dd,                                           \
            0,                                                          \
            sizeof (_ddenominator_dd));                                 \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        float _d = _logp[_logpri] - _max;                               \
        float _expd = fastexp (_d);                                     \
        _denominator += _expd;                                          \
        _ddenominator_dlogbeta += _expd * _d;                           \
                                                                        \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _ddenominator_dloga[_logprn][_logpri] -=                    \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _invrank *                  \
                                            _c (_logprn) *              \
                                            _d (_logprn, _logpri);      \
            _ddenominator_dloga[_logprn][_maxindex] +=                  \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _invrank *                  \
                                            _c (_logprn) *              \
                                            _d (_logprn, _maxindex);    \
            _ddenominator_dc[_logprn][_logpri] -=                       \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _invrank *                  \
                                            _d (_logprn, _logpri);      \
            _ddenominator_dc[_logprn][_maxindex] +=                     \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _invrank *                  \
                                            _d (_logprn, _maxindex);    \
            _ddenominator_dd[_logprn][_logpri] -=                       \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _invrank *                  \
                                            _c (_logprn);               \
            _ddenominator_dd[_logprn][_maxindex] +=                     \
                                            _expd *                     \
                                            _a[_logprn] *               \
                                            _beta *                     \
                                            _invrank *                  \
                                            _c (_logprn);               \
          }                                                             \
      }                                                                 \
                                                                        \
    for (unsigned int _logpri = 0; _logpri < _n_labels; ++_logpri)      \
      {                                                                 \
        for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)      \
          {                                                             \
            _dloga (_logprn) -= _pzj *                                  \
              _ddenominator_dloga[_logprn][_logpri] / _denominator;     \
            _dc (_logprn) -= _pzj *                                     \
              _ddenominator_dc[_logprn][_logpri] / _denominator;        \
            _dd (_logprn, _logpri) += _pzj *                            \
              -_ddenominator_dd[_logprn][_logpri] / _denominator;       \
          }                                                             \
      }                                                                 \
                                                                        \
    for (unsigned int _logprn = 0; _logprn < _rank; ++_logprn)          \
      {                                                                 \
        _dloga (_logprn) -= _pzj *                                      \
                           _a[_logprn] *                                \
                           _beta *                                      \
                           _invrank *                                   \
                           _c (_logprn) *                               \
                           _d (_logprn, _rating);                       \
        _dloga (_logprn) += _pzj *                                      \
                            _a[_logprn] *                               \
                            _beta *                                     \
                            _invrank *                                  \
                            _c (_logprn) *                              \
                            _d (_logprn, _maxindex);                    \
        _dc (_logprn) -= _pzj *                                         \
                         _a[_logprn] *                                  \
                         _beta *                                        \
                         _invrank *                                     \
                         _d (_logprn, _rating);                         \
        _dc (_logprn) += _pzj *                                         \
                         _a[_logprn] *                                  \
                         _beta *                                        \
                         _invrank *                                     \
                         _d (_logprn, _maxindex);                       \
        _dd (_logprn, _rating) -= _pzj *                                \
                                  _a[_logprn] *                         \
                                  _beta *                               \
                                  _invrank *                            \
                                  _c (_logprn);                         \
        _dd (_logprn, _maxindex) += _pzj *                              \
                                    _a[_logprn] *                       \
                                    _beta *                             \
                                    _invrank *                          \
                                    _c (_logprn);                       \
      }                                                                 \
  })
